Load Libraries

Description: Load all required packages for data retrieval, manipulation, and plotting.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(httr)
library(jsonlite)
## 
## Attaching package: 'jsonlite'
## 
## The following object is masked from 'package:purrr':
## 
##     flatten
library(ggplot2)
library(purrr)
library(zoo)      
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(plotly)   
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:httr':
## 
##     config
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

Download USGS earthquake data for Trident region

Description: Define the Trident bounding box and time window, then query the USGS FDSN API for all earthquakes in that region and time period.

# 1. Define Trident region and time window
trident_bbox <- list(
  minlatitude  = 58.0,
  maxlatitude  = 58.5,
  minlongitude = -155.4,
  maxlongitude = -154.8
)

starttime <- "2010-01-01"
endtime   <- "2025-11-22"   # update as needed

# 2. Build and send request to USGS
base_url <- "https://earthquake.usgs.gov/fdsnws/event/1/query"

resp <- GET(
  url = base_url,
  query = c(
    format = "geojson",
    starttime = starttime,
    endtime   = endtime,
    trident_bbox,
    eventtype = "earthquake",
    minmagnitude = 0
  )
)

stop_for_status(resp)

raw_json <- content(resp, "text", encoding = "UTF-8")
eq <- fromJSON(raw_json, flatten = TRUE)

Build event-level earthquake table

Description: Convert the raw GeoJSON features into a tidy event-level table with time, magnitude, location, and depth (in km and miles).

quakes <- eq$features %>%
  as_tibble() %>%
  transmute(
    time = as.POSIXct(properties.time / 1000,
                      origin = "1970-01-01", tz = "UTC"),
    mag   = properties.mag,
    place = properties.place,
    lon   = purrr::map_dbl(geometry.coordinates, 1),
    lat   = purrr::map_dbl(geometry.coordinates, 2),
    depth_km    = purrr::map_dbl(geometry.coordinates, 3),      # depth from USGS
    depth_miles = depth_km * 0.621371                          # depth for public-facing labels
  ) %>%
  arrange(time)

glimpse(quakes)
## Rows: 10,408
## Columns: 7
## $ time        <dttm> 2010-01-01 11:28:23, 2010-01-02 07:34:33, 2010-01-02 18:4…
## $ mag         <dbl> 1.0, 1.0, 1.2, 1.6, 1.2, 1.7, 1.4, 1.6, 2.2, 1.2, 1.0, 1.6…
## $ place       <chr> "99 km NNW of Karluk, Alaska", "86 km NW of Karluk, Alaska…
## $ lon         <dbl> -154.9000, -155.3670, -154.9041, -155.3792, -154.9891, -15…
## $ lat         <dbl> 58.4336, 58.1697, 58.4343, 58.1675, 58.2845, 58.1774, 58.1…
## $ depth_km    <dbl> 1.0, 5.3, 1.0, 1.0, 2.7, 4.6, 2.8, 6.5, 1.0, 3.0, 3.2, 1.0…
## $ depth_miles <dbl> 0.6213710, 3.2932663, 0.6213710, 0.6213710, 1.6777017, 2.8…

Aggregate to daily seismic metrics

Description: Summarize event-level data into daily metrics that will feed into the unrest index (daily event count, magnitude distribution, and shallow-earthquake percentage). Negative depths are set to 0 km, and “shallow” is defined as < 3 km.

quakes_daily <- quakes %>%
  mutate(
    depth_km = ifelse(depth_km < 0, 0, depth_km),   # clamp negative depths to 0 km
    date     = as.Date(time)
  ) %>%
  group_by(date) %>%
  summarise(
    n_events    = n(),                                # daily number of earthquakes
    median_mag  = median(mag, na.rm = TRUE),          # typical quake size that day
    p90_mag     = quantile(mag, 0.9, na.rm = TRUE),   # larger quake size indicator
    shallow_pct = mean(depth_km <= 3, na.rm = TRUE),  # % of shallow quakes (< 3 km)
    .groups = "drop"                                  # return ungrouped dataframe
  )

head(quakes_daily)
## # A tibble: 6 × 5
##   date       n_events median_mag p90_mag shallow_pct
##   <date>        <int>      <dbl>   <dbl>       <dbl>
## 1 2010-01-01        1        1      1          1    
## 2 2010-01-02        3        1.2    1.52       0.667
## 3 2010-01-09        1        1.2    1.2        1    
## 4 2010-01-12        1        1.7    1.7        0    
## 5 2010-01-14        1        1.4    1.4        1    
## 6 2010-01-15        2        1.9    2.14       0.5

Baseline configuration (fixed, rolling, anchored)

Description: Configure how the baseline is defined. • baseline_mode = “fixed” uses a pre-unrest historical window. • baseline_mode = “rolling” uses a rolling window (previous N years) for each date. • baseline_mode = “anchored” uses a 3/5/7-year or historical window before a chosen reference date t₀.

# -------------------------------------------------------------------
# Baseline configuration
# -------------------------------------------------------------------

# Choose baseline mode: "fixed", "rolling", or "anchored"
baseline_mode <- "anchored"

# For rolling baseline (not the main focus right now)
baseline_window_years <- 3   # previous N years per date if baseline_mode == "rolling"

# Fixed (historical) baseline: user-defined quiet period (pre-unrest)
baseline_fixed_start <- as.Date("2015-01-01")
baseline_fixed_end   <- as.Date("2022-08-24")  # start of documented unrest at Trident

# Anchored baseline: history before a chosen reference date t0
anchored_ref_date     <- as.Date("2022-08-24")  # t0; can be any date of interest
anchored_span_years   <- 5                      # 3, 5, or 7 years
anchored_historical   <- FALSE                  # TRUE = full history up to t0 (ignore span)

Rolling baseline function (for optional rolling mode)

Description: For rolling mode, compute per-date baseline statistics using the previous N years of daily seismic metrics (excluding the current day).

compute_rolling_baseline <- function(quakes_daily,
                                     window_years = 3,
                                     min_days = 30) {
  dates <- quakes_daily$date
  
  purrr::map_dfr(dates, function(d) {
    start <- d %m-% years(window_years)  # window_years before date d
    ref <- quakes_daily %>%
      filter(date >= start, date < d)    # strictly before current day
    
    if (nrow(ref) < min_days) {
      tibble(
        date   = d,
        n_mu   = NA_real_,
        n_sd   = NA_real_,
        p90_mu = NA_real_,
        p90_sd = NA_real_,
        shp_mu = NA_real_,
        shp_sd = NA_real_
      )
    } else {
      tibble(
        date   = d,
        n_mu   = mean(ref$n_events,    na.rm = TRUE),
        n_sd   = sd(ref$n_events,      na.rm = TRUE),
        p90_mu = mean(ref$p90_mag,     na.rm = TRUE),
        p90_sd = sd(ref$p90_mag,       na.rm = TRUE),
        shp_mu = mean(ref$shallow_pct, na.rm = TRUE),
        shp_sd = sd(ref$shallow_pct,   na.rm = TRUE)
      )
    }
  })
}

Compute unrest score using selected baseline mode

Description: • Standardize each daily metric relative to its baseline (fixed, rolling, or anchored) using z-scores. • Convert z-scores to intensities (0–1), average them, and rescale to 0–100. • Assign a traffic-light status (Green/Yellow/Orange/Red) based on the unrest score.

intensity_from_z <- function(z) {
  x <- z / 3           # 3 sd ~= "max" intensity
  x <- pmin(pmax(x, 0), 1)
  return(x)
}

eps <- 1e-6  # to avoid divide-by-zero

if (baseline_mode == "fixed") {
  # ----------------------------------------------
  # Fixed baseline: single historical quiet period
  # ----------------------------------------------
  baseline_stats <- quakes_daily %>%
    filter(date >= baseline_fixed_start,
           date <  baseline_fixed_end) %>%
    summarise(
      n_mu    = mean(n_events,    na.rm = TRUE),
      n_sd    = sd(n_events,      na.rm = TRUE),
      p90_mu  = mean(p90_mag,     na.rm = TRUE),
      p90_sd  = sd(p90_mag,       na.rm = TRUE),
      shp_mu  = mean(shallow_pct, na.rm = TRUE),
      shp_sd  = sd(shallow_pct,   na.rm = TRUE)
    )
  
  quakes_scored <- quakes_daily %>%
    mutate(
      z_n   = (n_events    - baseline_stats$n_mu)   / (baseline_stats$n_sd   + eps),
      z_p90 = (p90_mag     - baseline_stats$p90_mu) / (baseline_stats$p90_sd + eps),
      z_shp = (shallow_pct - baseline_stats$shp_mu) / (baseline_stats$shp_sd + eps),
      
      I_n   = intensity_from_z(z_n),
      I_p90 = intensity_from_z(z_p90),
      I_shp = intensity_from_z(z_shp),
      
      unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
      status = case_when(
        unrest_score < 20 ~ "Green",
        unrest_score < 40 ~ "Yellow",
        unrest_score < 70 ~ "Orange",
        TRUE              ~ "Red"
      )
    )
  
} else if (baseline_mode == "rolling") {
  # ----------------------------------------------
  # Rolling baseline: previous N years per date
  # ----------------------------------------------
  baseline_rolling <- compute_rolling_baseline(
    quakes_daily,
    window_years = baseline_window_years
  )
  
  quakes_scored <- quakes_daily %>%
    left_join(baseline_rolling, by = "date") %>%
    mutate(
      z_n   = (n_events    - n_mu)   / (n_sd   + eps),
      z_p90 = (p90_mag     - p90_mu) / (p90_sd + eps),
      z_shp = (shallow_pct - shp_mu) / (shp_sd + eps),
      
      I_n   = intensity_from_z(z_n),
      I_p90 = intensity_from_z(z_p90),
      I_shp = intensity_from_z(z_shp),
      
      unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
      status = case_when(
        unrest_score < 20 ~ "Green",
        unrest_score < 40 ~ "Yellow",
        unrest_score < 70 ~ "Orange",
        TRUE              ~ "Red"
      )
    )
  
} else if (baseline_mode == "anchored") {
  # ----------------------------------------------
  # Anchored baseline: history before a reference date t0
  # ----------------------------------------------
  # Determine baseline window relative to anchored_ref_date
  if (anchored_historical) {
    baseline_start <- min(quakes_daily$date, na.rm = TRUE)
  } else {
    baseline_start <- anchored_ref_date %m-% years(anchored_span_years)
  }
  baseline_end <- anchored_ref_date
  
  baseline_stats <- quakes_daily %>%
    filter(date >= baseline_start,
           date <  baseline_end) %>%
    summarise(
      n_mu    = mean(n_events,    na.rm = TRUE),
      n_sd    = sd(n_events,      na.rm = TRUE),
      p90_mu  = mean(p90_mag,     na.rm = TRUE),
      p90_sd  = sd(p90_mag,       na.rm = TRUE),
      shp_mu  = mean(shallow_pct, na.rm = TRUE),
      shp_sd  = sd(shallow_pct,   na.rm = TRUE)
    )
  
  quakes_scored <- quakes_daily %>%
    mutate(
      z_n   = (n_events    - baseline_stats$n_mu)   / (baseline_stats$n_sd   + eps),
      z_p90 = (p90_mag     - baseline_stats$p90_mu) / (baseline_stats$p90_sd + eps),
      z_shp = (shallow_pct - baseline_stats$shp_mu) / (baseline_stats$shp_sd + eps),
      
      I_n   = intensity_from_z(z_n),
      I_p90 = intensity_from_z(z_p90),
      I_shp = intensity_from_z(z_shp),
      
      unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
      status = case_when(
        unrest_score < 20 ~ "Green",
        unrest_score < 40 ~ "Yellow",
        unrest_score < 70 ~ "Orange",
        TRUE              ~ "Red"
      )
    )
}

Smooth unrest score with a 7-day moving average

Description: To reduce noise from day-to-day fluctuations in seismicity, a 7-day rolling mean is applied to the unrest index. This smoothing highlights broader trends in volcanic activity, making the seismic state easier to interpret in dashboards and public-facing visualizations. The raw daily score is preserved for internal analysis, while the averaged score (unrest_score_7d) is used for plotting.

quakes_scored <- quakes_scored %>%
  arrange(date) %>%
  mutate(
    unrest_score_7d = as.numeric(zoo::rollapply(
      unrest_score,
      width = 7,
      FUN = mean,
      align = "right",
      fill = NA,
      na.rm = TRUE
    ))
  )

summary(quakes_scored$unrest_score_7d)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.9969 10.3291 13.8954 16.9065 20.1364 50.8806       6

Baseline mean for plotting (anchored to t₀)

Description: Compute the mean of the 7-day smoothed unrest score over the chosen baseline window (3/5/7 years or historical before the reference date). This will be shown as a dashed horizontal line on the plot.

if (baseline_mode == "fixed") {
  baseline_start_plot <- baseline_fixed_start
  baseline_end_plot   <- baseline_fixed_end
  
} else if (baseline_mode == "anchored") {
  if (anchored_historical) {
    baseline_start_plot <- min(quakes_scored$date, na.rm = TRUE)
  } else {
    baseline_start_plot <- anchored_ref_date %m-% years(anchored_span_years)
  }
  baseline_end_plot <- anchored_ref_date
  
} else {
  baseline_start_plot <- NA
  baseline_end_plot   <- NA
}

if (!is.na(baseline_start_plot)) {
  baseline_mean_score <- quakes_scored %>%
    filter(date >= baseline_start_plot,
           date <  baseline_end_plot) %>%
    summarise(mu = mean(unrest_score_7d, na.rm = TRUE)) %>%
    pull(mu)
} else {
  baseline_mean_score <- NA_real_
}

baseline_start_plot
## [1] "2017-08-24"
baseline_end_plot
## [1] "2022-08-24"
baseline_mean_score
## [1] 12.90601

Plot the Trident Seismic Unrest Score

Description: Visualize the 7-day smoothed unrest score over time with color-coded status, a vertical line marking the anchored reference date, and a dashed horizontal line showing the baseline mean over the chosen window.

subtitle_text <- dplyr::case_when(
  baseline_mode == "fixed" ~ paste0(
    "Fixed baseline: ", baseline_fixed_start, " to ", baseline_fixed_end, " (pre-unrest)"
  ),
  baseline_mode == "rolling" ~ paste0(
    "Rolling baseline: previous ", baseline_window_years, " years per day"
  ),
  baseline_mode == "anchored" & anchored_historical ~ paste0(
    "Anchored historical baseline up to ", anchored_ref_date
  ),
  baseline_mode == "anchored" ~ paste0(
    "Anchored baseline: ", anchored_span_years, "-year window before ", anchored_ref_date
  ),
  TRUE ~ ""
)

# Vertical line date: baseline_fixed_end for fixed, anchored_ref_date for anchored
vline_date <- if (baseline_mode == "fixed") {
  baseline_fixed_end
} else if (baseline_mode == "anchored") {
  anchored_ref_date
} else {
  NA
}

p <- ggplot(
  quakes_scored,
  aes(
    x = date,
    y = unrest_score_7d,
    color = status,
    text = paste0(
      "Date: ", date,
      "<br>Score (7d): ", round(unrest_score_7d, 1),
      "<br>Daily events: ", n_events,
      "<br>p90 mag: ", round(p90_mag, 2),
      "<br>Shallow %: ", round(shallow_pct * 100, 1), "%"
    )
  )
) +
  geom_line(na.rm = TRUE) +
  geom_vline(
    xintercept = as.numeric(vline_date),
    linetype   = "dashed"
  ) +
  {
    if (!is.na(baseline_mean_score)) {
      geom_hline(yintercept = baseline_mean_score,
                 linetype = "dashed")
    } else {
      NULL
    }
  } +
# Baseline label annotation (rounded to 2 decimal places, placed above the line)
{
  if (!is.na(baseline_mean_score)) {
    annotate(
      "text",
      x = max(quakes_scored$date, na.rm = TRUE),
      y = baseline_mean_score + 2,  # 2 units above the line on the 0–100 scale
      label = paste0("Baseline \u2248 ", round(baseline_mean_score, 2)),
      hjust = -0.1,
      vjust = 0.5,
      size = 3.5,
      color = "black"
    )
  } else {
    NULL
  }
} +
  labs(
    title = "Trident Volcano – 7-Day Seismic Unrest Score",
    subtitle = subtitle_text,
    x = "Date",
    y = "Unrest Score (0–100, 7-day mean)",
    color = "Status"
  ) +
  scale_color_manual(
    values = c(
      "Green"  = "#00CC00",
      "Yellow" = "#FFD700",
      "Orange" = "#FF9900",
      "Red"    = "#CC0000"
    ),
    breaks = c("Red", "Orange", "Yellow", "Green")
  ) +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "#F0F0F0", color = NA),
    plot.background  = element_rect(fill = "#F0F0F0", color = NA),
    panel.grid.major = element_line(color = "grey80"),
    panel.grid.minor = element_line(color = "grey90"),
    clip = "off"
  )

Interactive hover-enabled unrest visualization

Description: The final plot is rendered as an interactive HTML widget using plotly::ggplotly(). Hovering over the line (and points) reveals detailed values for each date, including the 7-day unrest score, daily earthquake count, magnitude distribution, and shallow earthquake percentage.

# Add points to make individual days easy to see in the interactive view
p <- p + geom_point(size = 0.5, alpha = 0.9)

ggplotly(p, tooltip = "text")
## Warning in scale_x_date(): A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.
## Warning in plot_theme(list(theme = theme), default = default): The `clip` theme
## element is not defined in the element hierarchy.